stanEx1
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}stan basic
normal distribution
ex1-0.stan
normal distribution
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
real y1;
y1=normal_rng(m,s);
}mdl=cmdstan_model('./ex1-0.stan')
y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.66 1.74 2.20 2.17 2.40 3.42par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y') 
data=list(N=length(y),y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -32.4689
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 9 -2.00847 0.000241842 0.000187564 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.3 seconds.mle## variable estimate
## lp__ -2.01
## m 2.17
## s 0.74
## y1 1.22mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -3.35 -3.01 1.04 0.75 -5.46 -2.34 1.00 803 925
## m 2.16 2.16 0.30 0.28 1.69 2.62 1.00 1091 918
## s 0.91 0.86 0.25 0.21 0.60 1.37 1.00 1380 1165
## y1 2.14 2.15 0.99 0.91 0.55 3.72 1.00 1907 1879mcmc$metadata()$stan_variables## [1] "lp__" "m" "s" "y1"mcmc$metadata()$model_params## [1] "lp__" "m" "s" "y1"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -3.35 -3.01 1.04 0.75 -5.46 -2.34 1.00 803 925
## m 2.16 2.16 0.30 0.28 1.69 2.62 1.00 1091 918
## s 0.91 0.86 0.25 0.21 0.60 1.37 1.00 1380 1165
## y1 2.14 2.15 0.99 0.91 0.55 3.72 1.00 1907 1879
poisson distribution
y=rpois(20,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.75 2.50 2.95 3.25 10.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)ex2-1.stan
in case fit poisson dist. to normal dist.
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
real y1;
y1=normal_rng(m,s);
}# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -33.8038
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -26.7663 0.000851482 0.00122998 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -26.77
## m 2.95
## s 2.31
## y1 8.05mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -25.88 -25.59 1.02 0.77 -27.97 -24.88 1.00 786 807
## m 2.96 2.94 0.58 0.58 2.04 3.91 1.00 1065 870
## s 2.53 2.48 0.43 0.40 1.93 3.31 1.00 1610 1257
## y1 2.96 2.88 2.63 2.52 -1.43 7.37 1.00 1932 1974
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-2.stan
fit to poisson dist.
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
int y1;
y1=poisson_rng(l);
}mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -20.7726
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 4.82651 0.00029093 0.000121597 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 4.83
## l 2.95
## y1 4.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 5.44 5.71 0.66 0.28 4.05 5.92 1.00 833 1284
## l 3.02 3.00 0.38 0.37 2.39 3.67 1.00 819 938
## y1 3.00 3.00 1.79 1.48 0.00 6.00 1.00 1710 1869
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
difference test
ex3.stan
mean difference of 2 normal distributions
data{
int N;
vector[N] a;
vector[N] b;
}
parameters{
real ma;
real<lower=0> sa;
real mb;
real<lower=0> sb;
}
model{
a~normal(ma,sa);
b~normal(mb,sb);
}
generated quantities{
real d;
d=mb-ma;
}n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)
mdl=cmdstan_model('./ex3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -33482.6
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -33.3995 0.000938568 0.00156273 0.9537 0.9537 59
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -33.40
## ma 10.14
## sa 0.78
## mb 11.13
## sb 2.49
## d 0.99mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -34.75 -34.42 1.42 1.26 -37.52 -33.03 1.00 799 1160
## ma 10.15 10.15 0.19 0.20 9.83 10.45 1.00 1892 1204
## sa 0.87 0.84 0.16 0.14 0.65 1.17 1.00 2124 1341
## mb 11.12 11.13 0.59 0.56 10.10 12.08 1.01 476 619
## sb 2.74 2.67 0.49 0.46 2.08 3.67 1.00 2226 1544
## d 0.97 0.98 0.62 0.60 -0.06 1.98 1.00 537 709
d=mcmc$draws('d')
d## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
##
## chain
## iteration 1 2 3 4
## 1 2.215 0.92 1.68 0.76
## 2 0.068 0.55 2.43 0.98
## 3 0.659 0.80 1.63 1.19
## 4 1.960 0.47 1.51 2.31
## 5 1.953 0.39 -0.41 2.17
##
## # ... with 495 more iterationsmean(d>0)## [1] 0.935normal regression
ex4-1.stan
single normal regression
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -50125.8
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 34 -58.8763 0.00241247 0.000958411 0.945 0.945 90
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -58.88
## b0 15.00
## b1 2.95
## s 11.52mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -58.16 -57.66 1.55 1.09 -61.15 -56.62 1.01 302 225
## b0 15.18 15.04 6.50 5.60 5.25 26.29 1.02 176 138
## b1 2.95 2.95 0.10 0.09 2.78 3.11 1.01 229 241
## s 13.22 12.81 2.47 2.16 9.97 17.75 1.01 746 594
b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)
ex4-2.stan
single normal regression, prediction from new explanatory variable x1
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N1] m;
vector[N1] y1;
for(i in 1:N1){
m[i]=b0+b1*x1[i];
y1[i]=normal_rng(m[i],s);
}
}# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex4-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -27844.4
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 27 -58.8763 0.0076696 0.00260617 1 1 71
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -58.88
## b0 14.99
## b1 2.95
## s 11.52
## m[1] 22.35
## m[2] 52.97
## m[3] 83.60
## m[4] 114.23
## m[5] 144.85
## m[6] 175.48
## m[7] 206.11
## m[8] 236.74
## m[9] 267.36
## m[10] 297.99
## y1[1] 13.58
## y1[2] 49.62
## y1[3] 88.66
## y1[4] 124.27
## y1[5] 151.62
## y1[6] 157.93
## y1[7] 206.78
## y1[8] 227.16
## y1[9] 303.72
## y1[10] 303.28mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -58.04 -57.76 1.29 1.09 -60.57 -56.62 1.01 445 595
## b0 14.69 14.65 5.61 5.40 5.48 23.77 1.02 249 180
## b1 2.95 2.95 0.09 0.09 2.80 3.10 1.01 326 337
## s 13.07 12.79 2.40 2.21 9.93 17.33 1.00 1100 881
## m[1] 22.05 22.04 5.41 5.22 13.20 30.83 1.02 250 183
## m[2] 52.72 52.68 4.64 4.38 45.01 60.25 1.02 253 185
## m[3] 83.40 83.34 3.95 3.80 76.81 89.93 1.02 268 210
## m[4] 114.07 113.99 3.39 3.28 108.48 119.63 1.01 313 323
## m[5] 144.74 144.68 3.03 2.90 139.85 149.74 1.00 450 580
## m[6] 175.42 175.39 2.95 2.84 170.64 180.24 1.00 967 924
## m[7] 206.09 206.07 3.17 3.08 200.69 211.20 1.00 2319 1649
## m[8] 236.76 236.72 3.64 3.38 230.58 242.73 1.01 2458 1416
## m[9] 267.43 267.36 4.27 3.95 260.30 274.49 1.01 1869 1405
## m[10] 298.11 298.05 5.01 4.82 289.95 306.33 1.00 1221 1348
## y1[1] 21.81 21.85 14.24 13.78 -1.53 44.47 1.00 979 1344
## y1[2] 52.40 52.33 14.02 13.50 29.62 75.10 1.00 1010 1136
## y1[3] 83.15 83.12 13.82 13.27 60.58 104.87 1.00 1270 1614
## y1[4] 114.63 114.20 13.73 13.61 93.29 137.27 1.00 1617 1858
## y1[5] 144.31 144.88 13.76 13.42 121.85 166.26 1.00 1833 1720
## y1[6] 175.56 175.58 13.52 12.70 153.84 197.48 1.00 1873 2031
## y1[7] 206.43 206.31 13.81 12.91 183.85 228.30 1.00 1800 1984
## y1[8] 236.46 236.72 13.80 13.16 213.56 259.08 1.00 1823 1855
## y1[9] 266.85 266.52 13.97 13.34 244.01 290.08 1.00 1823 1741
## y1[10] 298.22 298.30 14.38 13.73 274.54 321.91 1.00 2153 1598
m=mcmc$draws('m')
summary(m)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 m[1] 22.1 22.0 5.41 5.22 13.2 30.8 1.02 250. 184.
## 2 m[2] 52.7 52.7 4.64 4.38 45.0 60.3 1.02 254. 186.
## 3 m[3] 83.4 83.3 3.95 3.80 76.8 89.9 1.02 269. 211.
## 4 m[4] 114. 114. 3.39 3.28 108. 120. 1.01 313. 324.
## 5 m[5] 145. 145. 3.03 2.90 140. 150. 1.00 451. 580.
## 6 m[6] 175. 175. 2.95 2.84 171. 180. 1.00 968. 924.
## 7 m[7] 206. 206. 3.17 3.08 201. 211. 1.00 2319. 1649.
## 8 m[8] 237. 237. 3.64 3.38 231. 243. 1.01 2459. 1417.
## 9 m[9] 267. 267. 4.27 3.95 260. 274. 1.01 1870. 1406.
## 10 m[10] 298. 298. 5.01 4.82 290. 306. 1.00 1221. 1349.smm=summary(m)
y1=mcmc$draws('y1')
summary(y1)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 y1[1] 21.8 21.8 14.2 13.8 -1.53 44.5 1.00 979. 1344.
## 2 y1[2] 52.4 52.3 14.0 13.5 29.6 75.1 1.00 1011. 1137.
## 3 y1[3] 83.2 83.1 13.8 13.3 60.6 105. 1.00 1270. 1615.
## 4 y1[4] 115. 114. 13.7 13.6 93.3 137. 0.999 1617. 1859.
## 5 y1[5] 144. 145. 13.8 13.4 122. 166. 1.00 1833. 1720.
## 6 y1[6] 176. 176. 13.5 12.7 154. 197. 1.00 1873. 2032.
## 7 y1[7] 206. 206. 13.8 12.9 184. 228. 1.00 1800. 1985.
## 8 y1[8] 236. 237. 13.8 13.2 214. 259. 1.00 1823. 1855.
## 9 y1[9] 267. 267. 14.0 13.3 244. 290. 1.00 1823. 1741.
## 10 y1[10] 298. 298. 14.4 13.7 275. 322. 1.00 2153. 1599.smy=summary(y1)
xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)
par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=m),xy,col='black')+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')
ex4-3.stan
single normal regression, prediction from original explanatory variable x
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m;
vector[N] y1;
for(i in 1:N){
m[i]=b0+b1*x[i];
y1[i]=normal_rng(m[i],s);
}
}data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-3.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -58.15 -57.84 1.38 1.19 -60.76 -56.62 1.01 281 272
## b0 14.44 15.07 6.47 6.43 3.59 24.17 1.01 167 143
## b1 2.96 2.95 0.10 0.10 2.80 3.13 1.01 219 216
## s 13.13 12.83 2.40 2.27 9.75 17.57 1.01 808 841
## m[1] 146.73 146.95 3.18 3.11 141.35 151.75 1.01 279 258
## m[2] 265.58 265.60 4.12 3.90 258.85 272.30 1.00 1452 1279
## m[3] 200.49 200.59 3.01 2.75 195.28 205.48 1.00 1535 901
## m[4] 215.81 215.85 3.16 2.93 210.44 221.07 1.00 2344 1175
## m[5] 26.49 27.10 6.10 6.07 16.33 35.67 1.01 167 136
## m[6] 278.42 278.53 4.44 4.21 271.27 285.67 1.00 1169 1150
## m[7] 298.14 298.27 4.97 4.78 289.99 306.24 1.00 907 1082
## m[8] 151.22 151.43 3.12 3.03 145.95 156.18 1.01 301 275
## m[9] 44.95 45.51 5.55 5.54 35.59 53.21 1.01 167 133
## m[10] 150.42 150.63 3.13 3.04 145.12 155.40 1.01 296 275
## m[11] 172.14 172.24 2.96 2.84 167.12 176.99 1.00 511 507
## m[12] 272.57 272.62 4.29 4.07 265.63 279.59 1.00 1273 1212
## m[13] 101.53 101.88 4.04 3.96 94.47 107.61 1.01 186 160
## m[14] 131.78 132.06 3.41 3.28 125.96 137.03 1.01 229 223
## m[15] 281.29 281.42 4.51 4.29 273.98 288.69 1.00 1123 1065
## m[16] 84.90 85.34 4.45 4.35 77.24 91.54 1.01 177 144
## m[17] 282.64 282.78 4.55 4.34 275.27 290.10 1.00 1102 1065
## m[18] 296.20 296.35 4.91 4.72 288.15 304.24 1.00 928 1073
## m[19] 21.81 22.44 6.24 6.21 11.40 31.25 1.01 167 139
## m[20] 24.46 25.09 6.16 6.11 14.17 33.76 1.01 167 138
## y1[1] 146.81 147.07 13.85 13.23 123.78 169.41 1.01 1842 1850
## y1[2] 265.51 265.72 13.71 13.29 242.70 288.06 1.00 1979 1794
## y1[3] 200.03 200.42 13.99 13.05 177.31 222.05 1.00 1879 1825
## y1[4] 216.19 216.34 13.97 13.53 193.34 238.59 1.00 1912 1848
## y1[5] 26.84 27.16 14.87 14.42 2.60 50.41 1.00 1094 1540
## y1[6] 279.31 278.90 13.97 13.30 257.20 302.44 1.00 1848 1988
## y1[7] 298.62 298.46 14.27 13.26 275.75 321.73 1.00 2039 1635
## y1[8] 151.51 151.78 13.79 13.20 128.37 173.46 1.00 1876 1761
## y1[9] 45.56 45.56 14.57 14.60 21.29 68.85 1.00 751 1492
## y1[10] 150.44 150.36 13.22 12.44 128.80 172.38 1.00 2018 1773
## y1[11] 171.94 172.02 13.37 12.94 149.26 193.84 1.00 2187 1957
## y1[12] 272.22 271.98 14.35 13.70 248.83 295.12 1.00 1955 1895
## y1[13] 101.34 101.77 13.72 13.06 78.58 122.87 1.00 1590 1813
## y1[14] 132.00 132.33 13.88 13.23 108.43 154.81 1.00 2008 1820
## y1[15] 281.55 281.64 13.78 12.87 258.73 303.51 1.00 2025 1673
## y1[16] 84.90 85.07 14.42 13.72 61.50 107.17 1.00 1192 1370
## y1[17] 283.34 283.04 14.07 12.95 260.59 307.26 1.00 1861 1897
## y1[18] 295.94 295.57 14.26 14.04 272.46 319.17 1.00 1768 1933
## y1[19] 21.79 21.80 14.73 14.31 -3.10 45.41 1.00 863 1135
## y1[20] 24.35 24.63 14.85 14.23 0.00 47.42 1.00 645 895
y1=mcmc$draws('y1')
smy=summary(y1)
par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
estimate correlation
ex4-41.stan
use covariance matrix and multinormal distribution
data{
int N;
array[N] vector[2] y;
}
parameters{
vector[2] m;
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
y~multi_normal(m,cv);
}ex4-42.stan
use covariance matrix and wishart distribution
data{
int N;
matrix[2,2] dp;
}
parameters{
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
dp~wishart(N-1,cv);
}ex4-43.stan
use correlation matrix and wishart distribution
data{
int N;
int M;
matrix[M,M] dp;
}
parameters{
vector<lower=0>[M] s;
corr_matrix[M] r;
}
transformed parameters{
cov_matrix[M] cv;
cv=quad_form_diag(r,s);
}
model{
dp~wishart(N-1,cv);
}ex4-44.stan use covariance matrix and multinormal distribution
data{
int N;
int K;
array[N] vector[K] y;
}
parameters{
vector[K] m;
cov_matrix[K] cov;
}
model{
y~multi_normal(m,cov);
}library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)## [,1] [,2]
## [1,] 1.213 0.674
## [2,] 0.674 0.902cor(y)## [,1] [,2]
## [1,] 1.000 0.644
## [2,] 0.644 1.000plot(y)
#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -132.14
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -14.5172 5.04246e-05 0.00107316 0.3566 1 28
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.52
## m[1] -0.02
## m[2] 0.24
## s[1] 1.07
## s[2] 0.93
## r 0.64
## cv[1,1] 1.15
## cv[2,1] 0.64
## cv[1,2] 0.64
## cv[2,2] 0.86mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -18.41 -18.02 1.73 1.53 -21.71 -16.35 1.00 812 1202
## m[1] -0.03 -0.04 0.26 0.25 -0.45 0.43 1.01 1054 1071
## m[2] 0.24 0.24 0.24 0.22 -0.14 0.62 1.00 1276 1050
## s[1] 1.18 1.16 0.21 0.19 0.89 1.56 1.00 1305 1236
## s[2] 1.03 1.00 0.19 0.17 0.78 1.36 1.00 1333 1505
## r 0.58 0.61 0.16 0.15 0.26 0.80 1.01 427 660
## cv[1,1] 1.44 1.34 0.55 0.44 0.80 2.45 1.00 1305 1236
## cv[2,1] 0.75 0.68 0.39 0.32 0.25 1.42 1.01 517 791
## cv[1,2] 0.75 0.68 0.39 0.32 0.25 1.42 1.01 517 791
## cv[2,2] 1.10 1.01 0.44 0.35 0.61 1.85 1.00 1333 1505
#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n
mdl=cmdstan_model('./ex4-42.stan')
data=list(N=n,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -26.1847
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -14.7659 0.000218707 1.49717e-05 1 1 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.77
## s[1] 1.10
## s[2] 0.95
## r 0.64
## cv[1,1] 1.21
## cv[2,1] 0.67
## cv[1,2] 0.67
## cv[2,2] 0.90mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.44 -17.14 1.22 1.03 -19.81 -16.11 1.00 755 1168
## s[1] 1.19 1.16 0.21 0.19 0.91 1.59 1.00 1223 1318
## s[2] 1.03 1.00 0.19 0.17 0.77 1.37 1.00 1221 1190
## r 0.60 0.61 0.15 0.14 0.33 0.82 1.01 532 721
## cv[1,1] 1.47 1.35 0.56 0.44 0.83 2.52 1.00 1223 1318
## cv[2,1] 0.77 0.69 0.39 0.31 0.30 1.53 1.01 596 725
## cv[1,2] 0.77 0.69 0.39 0.31 0.30 1.53 1.01 596 725
## cv[2,2] 1.09 1.00 0.42 0.35 0.59 1.89 1.00 1221 1190
#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')
m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -99.3293
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -14.7659 0.000171539 0.000536133 1 1 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.77
## s[1] 1.10
## s[2] 0.95
## r[1,1] 1.00
## r[2,1] 0.64
## r[1,2] 0.64
## r[2,2] 1.00
## cv[1,1] 1.21
## cv[2,1] 0.67
## cv[1,2] 0.67
## cv[2,2] 0.90mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.70 -16.46 1.15 1.03 -18.90 -15.39 1.01 921 1403
## s[1] 1.18 1.16 0.20 0.19 0.90 1.56 1.00 1256 1360
## s[2] 1.02 0.99 0.17 0.16 0.78 1.32 1.01 1189 1214
## r[1,1] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## r[2,1] 0.59 0.61 0.15 0.15 0.31 0.81 1.00 889 1021
## r[1,2] 0.59 0.61 0.15 0.15 0.31 0.81 1.00 889 1021
## r[2,2] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## cv[1,1] 1.44 1.33 0.52 0.44 0.80 2.44 1.00 1256 1360
## cv[2,1] 0.74 0.68 0.36 0.30 0.29 1.41 1.00 831 764
## cv[1,2] 0.74 0.68 0.36 0.30 0.29 1.41 1.00 831 764
## cv[2,2] 1.06 0.98 0.38 0.31 0.60 1.76 1.01 1189 1214
mdl=cmdstan_model('./ex4-44.stan')
k=2
data=list(N=n,K=k,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -1569.73
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -14.5172 8.20179e-05 0.000961096 0.8987 0.8987 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.52
## m[1] -0.02
## m[2] 0.24
## cov[1,1] 1.15
## cov[2,1] 0.64
## cov[1,2] 0.64
## cov[2,2] 0.86mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.13 -15.79 1.67 1.49 -19.28 -14.09 1.00 815 1069
## m[1] -0.03 -0.02 0.29 0.27 -0.53 0.43 1.01 846 817
## m[2] 0.24 0.23 0.25 0.25 -0.18 0.65 1.01 1132 942
## cov[1,1] 1.78 1.62 0.78 0.60 0.92 3.14 1.00 1522 1339
## cov[2,1] 0.98 0.88 0.51 0.41 0.39 1.93 1.00 1460 1146
## cov[1,2] 0.98 0.88 0.51 0.41 0.39 1.93 1.00 1460 1146
## cov[2,2] 1.29 1.17 0.53 0.41 0.69 2.29 1.00 1406 1191
distributions
normal dist.
distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2ex1-0.stan
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
real y1;
y1=normal_rng(m,s);
}n=20
y=rnorm(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.029 0.835 1.536 1.473 2.105 2.867par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -38.7162
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -4.99263 0.000103972 9.41093e-05 0.9847 0.9847 18
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -4.99
## m 1.47
## s 0.78
## y1 2.39mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.22 -5.94 0.96 0.71 -8.18 -5.29 1.00 947 1127
## m 1.48 1.48 0.18 0.18 1.18 1.78 1.00 1363 1058
## s 0.84 0.83 0.14 0.14 0.64 1.10 1.00 1743 1208
## y1 1.46 1.45 0.89 0.84 0.03 2.93 1.00 1856 1888
Bernoulli dist.
The event with probablity p occur y=1, not occur y=0
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)ex1-1.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
y~bernoulli(p);
}
generated quantities{
real s;
s=(p*(1-p))^.5;
}n=10
y=sample(0:1,n,replace=T)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 1.0 1.0 0.8 1.0 1.0par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -5.0084
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 3 -5.00402 0.00127944 3.77503e-05 1 1 6
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -5.00
## p 0.80
## s 0.40mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.27 -7.00 0.69 0.35 -8.76 -6.75 1.00 909 859
## p 0.75 0.76 0.12 0.12 0.52 0.92 1.00 617 827
## s 0.41 0.43 0.07 0.07 0.27 0.50 1.00 566 827
geometric dist.
The event with probablity p occur after y trials(failure)
y~Ge(p), y[0,Inf), P(Y=y)=p(1−p)^y
E[y]=1/p, V[y]=(1-p)/p^2ex1-2.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
for (n in 1:N) {
target+=log(p)+y[n]*log(1-p);
}
}
generated quantities{
real m;
m=1/p;
real s;
s=((1-p)/p^2)^.5;
}n=20
y=rgeom(n,0.3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.75 1.50 3.05 6.00 10.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -45.641
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 3 -45.2724 0.00360832 0.000134258 0.9482 0.9482 5
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -45.27
## p 0.25
## m 4.05
## s 3.51mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -47.42 -47.17 0.65 0.30 -48.74 -46.95 1.00 903 1451
## p 0.25 0.25 0.05 0.05 0.18 0.33 1.00 747 952
## m 4.08 3.94 0.79 0.71 3.03 5.56 1.00 747 952
## s 3.55 3.40 0.80 0.72 2.48 5.03 1.00 747 952
negative binomial dist.
The event with probablity p occur a times after y-a trials
y~NB(a,p), y[0,Inf), P(Y=y)=p^a*(1−p)^(y-a)
E[y]=a*(1-p)/p, V[y]=a*(1-p)/p^2ex1-3.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
real<lower=1> a;
}
model{
y~neg_binomial(a,p);
}
generated quantities{
real m;
m=a*(1-p)/p;
real s;
s=(a*(1-p)/p^2)^.5;
int y1;
y1=neg_binomial_rng(a,p);
}n=50
a=3
y=rnbinom(n,a,0.3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 3.00 5.00 5.90 7.75 18.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -143.873
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -137.898 0.00382547 0.000802942 0.9311 0.9311 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -137.90
## p 0.41
## a 2.43
## m 3.47
## s 2.90
## y1 13.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -139.93 -139.57 1.11 0.76 -142.12 -138.86 1.00 481 541
## p 0.52 0.50 0.16 0.16 0.30 0.82 1.01 443 313
## a 3.03 2.89 0.89 0.88 1.81 4.69 1.01 448 455
## m 2.84 2.86 1.07 1.03 1.00 4.55 1.01 496 313
## s 2.44 2.39 0.85 0.82 1.11 3.84 1.01 466 326
## y1 5.78 5.00 4.34 4.45 1.00 14.00 1.00 1916 1799
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
binomial dist.
The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)ex2-3.stan
data{
int N;
int n;
array[N] int y;
}
parameters{
real<lower=0,upper=0> p;
}
model{
y~binomial(n,p);
}
generated quantities{
real m;
m=n*p;
real s;
s=(n*p*(1-p))^.5;
int y1;
y1=binomial_rng(n,p);
}n=20
n0=10
y=rbinom(n,n0,0.3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 2.0 3.0 3.1 4.0 6.0par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex2-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -144.737
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -123.82 0.00185335 0.000182482 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -123.82
## p 0.31
## m 3.10
## s 1.46
## y1 2.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -125.86 -125.60 0.69 0.33 -127.22 -125.36 1.00 1168 1356
## p 0.31 0.31 0.03 0.03 0.26 0.37 1.00 983 1175
## m 3.14 3.13 0.33 0.34 2.61 3.69 1.00 983 1175
## s 1.46 1.47 0.04 0.04 1.39 1.53 1.00 983 1175
## y1 3.11 3.00 1.51 1.48 1.00 6.00 1.00 1747 1964
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-31.stan
trial n is varied from each sample
data{
int N;
array[N] int n;
array[N] int y;
}
parameters{
real<lower=0> p;
}
model{
y~binomial(n,p);
}n=20
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex2-31.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -81.2821
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -66.326 0.00582635 0.00143967 1 1 6
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -66.33
## p 0.29mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -68.37 -68.10 0.68 0.27 -69.76 -67.90 1.00 1104 1270
## p 0.29 0.29 0.04 0.04 0.23 0.36 1.01 707 968
Poisson dist.
The event occur l times in unit range, the event occur y times.
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=lex1-4.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
real s;
s=l^.5;
int y1;
y1=poisson_rng(l);
}n=20
y=rpois(n,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.75 3.00 3.30 4.00 12.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-4.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -58.9394
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 12.7989 0.00515492 0.000519729 0.9927 0.9927 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 12.80
## l 3.30
## s 1.82
## y1 1.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 13.49 13.77 0.70 0.31 12.02 14.00 1.00 1172 1478
## l 3.37 3.35 0.41 0.41 2.72 4.09 1.01 881 1161
## s 1.83 1.83 0.11 0.11 1.65 2.02 1.01 881 1161
## y1 3.36 3.00 1.87 1.48 1.00 7.00 1.00 1881 1892
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
exponential dist.
The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2ex1-5.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> l;
}
model{
y~exponential(l);
}
generated quantities{
real m;
m=1/l;
real s;
s=1/l;
real y1;
y1=exponential_rng(l);
}n=20
y=rexp(n,2)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.035 0.166 0.443 0.473 0.666 1.337par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-5.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -12.4359
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -5.00923 0.000141713 6.78712e-05 0.9789 0.9789 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -5.01
## l 2.12
## m 0.47
## s 0.47
## y1 2.30mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -4.74 -4.47 0.70 0.32 -6.21 -4.24 1.00 1026 1339
## l 2.20 2.16 0.48 0.47 1.49 3.04 1.00 758 1014
## m 0.48 0.46 0.11 0.10 0.33 0.67 1.00 758 1014
## s 0.48 0.46 0.11 0.10 0.33 0.67 1.00 758 1014
## y1 0.47 0.30 0.50 0.33 0.03 1.44 1.00 1755 1759
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
gamma dist.
The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2ex1-6.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> a;
real<lower=0> l;
}
model{
y~gamma(a,l);
}n=20
y=rgamma(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.86 1.59 2.17 2.27 2.78 5.15par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-6.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -137.813
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -26.5906 0.000411502 0.000416302 0.92 0.92 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -26.59
## a 5.40
## l 2.38
## m 2.27
## s 0.98
## y1 1.42mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -24.84 -24.54 0.96 0.70 -26.84 -23.90 1.01 564 932
## a 6.12 5.92 1.73 1.66 3.60 9.36 1.01 400 402
## l 2.72 2.65 0.79 0.79 1.56 4.19 1.02 395 470
## m 2.26 2.25 0.21 0.20 1.95 2.63 1.00 1919 1555
## s 0.94 0.91 0.16 0.15 0.72 1.25 1.01 441 622
## y1 2.26 2.13 0.95 0.87 0.96 3.96 1.00 1786 1868
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
log normal dist.
logalithm of variable follows normal distribution
log y~N(m,s), y>0
E[y]=exp(m+s^2)
V[y]=exp(2*m+s^2)*(exp(s^2)-1)ex1-7.stan
data{
int N;
vector[N]<lower=0> y;
}
parameters{
real m0;
real<lower=0> s0;
}
model{
y~lognormal(m0,s0);
}n=30
y=rlnorm(n,0.5,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.39 0.95 2.05 3.26 4.67 14.38par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-7.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -77.4717
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 7 -41.985 8.65034e-05 0.000161723 0.9434 0.9434 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -41.98
## m0 0.73
## s0 0.98
## m 5.45
## s 4.28
## y1 7.25mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -43.03 -42.74 1.00 0.75 -45.07 -42.05 1.00 921 1444
## m0 0.74 0.75 0.19 0.19 0.42 1.04 1.00 966 971
## s0 1.04 1.03 0.15 0.14 0.83 1.30 1.00 1917 1437
## m 6.86 6.07 3.42 2.11 3.74 12.29 1.00 1250 1166
## s 5.68 4.86 3.35 2.01 2.69 11.04 1.00 1318 1175
## y1 3.47 1.98 4.90 1.77 0.38 12.15 1.00 2031 1958
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
beta dist.
use as prior of binomial dist parameter p
event with probabilty p0 occur a-1 times, do not occur b-1 times in a+b-2 trials
p~Be(a,b), p[0,1], p=p0^(a-1)*(1-p0)^(b-1)
E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)ex1-8.stan
data{
int N;
vector[N] p;
}
parameters{
real<lower=0> a;
real<lower=0> b;
}
model{
p~beta(a,b);
}n=20
p=rbeta(n,1,2)
summary(p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.007 0.114 0.264 0.333 0.446 0.961par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')
data=list(N=n,p=p)
mdl=cmdstan_model('./ex1-8.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -184.562
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 3.13888 7.26822e-05 7.60616e-07 1 1 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 3.14
## a 0.76
## b 1.42mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 2.40 2.69 0.95 0.68 0.43 3.30 1.00 906 1305
## a 0.86 0.84 0.22 0.21 0.53 1.25 1.01 701 1154
## b 1.66 1.61 0.47 0.46 0.97 2.50 1.01 632 877
Dirichlet dist.
use as prior of categorical dist parameter p
p~dir(p0), p[0,1]ex1-9.stan
data {
int N;
int K;
matrix[N, K] p;
}
parameters {
simplex[K] p0;
}
model {
for(i in 1:N){
p[i]~dirichlet(p0);
}
}library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)## V1 V2 V3
## Min. :0.001 Min. :0.000 Min. :0.000
## 1st Qu.:0.257 1st Qu.:0.028 1st Qu.:0.000
## Median :0.827 Median :0.066 Median :0.013
## Mean :0.622 Mean :0.258 Mean :0.120
## 3rd Qu.:0.944 3rd Qu.:0.456 3rd Qu.:0.105
## Max. :0.999 Max. :0.945 Max. :0.996boxplot(p)
data=list(N=n,K=ncol(p),p=p)
mdl=cmdstan_model('./ex1-9.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = 51.4479
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 94.4034 0.000122507 0.000508777 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 94.40
## p0[1] 0.56
## p0[2] 0.29
## p0[3] 0.15mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 89.69 90.01 1.06 0.72 87.68 90.64 1.01 928 1005
## p0[1] 0.55 0.55 0.06 0.05 0.45 0.64 1.00 1933 1392
## p0[2] 0.30 0.30 0.05 0.05 0.21 0.39 1.00 1676 1289
## p0[3] 0.15 0.15 0.03 0.03 0.11 0.21 1.00 1873 1494
